Foundations of Dissipative Particle Dynamics 
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We derive a mesoscopic modeling and simulation technique that is very close to the technique 
known as dissipative particle dynamics. The model is derived from molecular dynamics by means 
of a systematic coarse-graining procedure. This procedure links the forces between the dissipative 
particles to a hydrodynamic description of the underlying molecular dynamics (MD) particles. In 
particular, the dissipative particle forces are given directly in terms of the viscosity emergent from 
MD, while the interparticle energy transfer is similarly given by the heat conductivity derived from 
MD. In linking the microscopic and mesoscopic descriptions we thus rely on the macroscopic de- 
scription emergent from MD. Thus the rules governing our new form of dissipative particle dynamics 
reflect the underlying molecular dynamics; in particular all the underlying conservation laws carry 
over from the microscopic to the mesoscopic descriptions. We obtain the forces experienced by 
the dissipative particles together with an approximate form of the associated equilibrium distribu- 
tion. Whereas previously the dissipative particles were spheres of fixed size and mass, now they 
are defined as cells on a Voronoi lattice with variable masses and sizes. This Voronoi lattice arises 
naturally from the coarse-graining procedure which may be applied iteratively and thus represents 
a form of renormalisation-group mapping. It enables us to select any desired local scale for the 
mesoscopic description of a given problem. Indeed, the method may be used to deal with situations 
in which several different length scales are simultaneously present. We compare and contrast this 
new particulate model with existing continuum fluid dynamics techniques, which rely on a purely 
macroscopic and phenomenological approach. Simulations carried out with the present scheme show 
good agreement with theoretical predictions for the equilibrium behavior. 

Pacs numbers: 47.11. +j 47.10. +g 05.40.+j 



I. INTRODUCTION 



The non-equilibrium behavior of fluids continues to 
present a major challenge for both theory and numeri- 
cal simulation. In recent times, there has been growing 
interest in the study of so-called 'mesoscale' modeling 
and simulation methods, particularly for the description 
of the complex dynamical behavior of many kinds of soft 
condensed matter, whose properties have thwarted more 
conventional approaches. As an example, consider the 
case of complex fluids with many coexisting length and 
time scales, for which hydrodynamic descriptions are un- 
known and may not even exist. These kinds of fluids in- 
clude multi-phase flows, particulate and colloidal suspen- 
sions, polymers, and amphiphilic fluids, including emul- 
sions and microemulsions. Fluctuations and Brownian 
motion are often key features controlling their behavior. 

From the standpoint of traditional fluid dynamics, a 
general problem in describing such fluids is the lack of 
adequate continuum models. Such descriptions, which 
are usually based on simple conservation laws, approach 
the physical description from the macroscopic side, that 
is in a 'top down' manner, and have certainly proved 
successful for simple Newtonian fluids Q. For complex 
fluids, however, equivalent phenomenological representa- 



tions are usually unavailable and instead it is necessary 
to base the modeling approach on a microscopic (that is 
on a particulate) description of the system, thus work- 
ing from the bottom upwards, along the general lines of 
the program for statistical mechanics pioneered by Boltz- 
mann g] . Molecular dynamics (MD) presents itself as the 
most accurate and fundamental method || but it is far 
too computationally intensive to provide a practical op- 
tion for most hydrodynamic problems involving complex 
fluids. Over the last decade several alternative 'bottom 
up' strategies have therefore been introduced. Hydrody- 
namic lattice gases Q , which model the fluid as a discrete 
set of particles, represent a computationally efficient spa- 
tial and temporal discretization of the more conventional 
molecular dynamics. The lattice-Boltzmann method ||, 
originally derived from the lattice-gas paradigm by in- 
voking Boltzmann's Stosszahlansatz, represents an inter- 
mediate (fluctuationless) approach between the top-down 
(continuum) and bottom-up (particulate) strategies, in- 
sofar as the basic entity in such models is a single particle 
distribution function; but for interacting systems even 
these lattice-Boltzmann methods can be subdivided into 
bottom- up H] and top-down models {?J . 

A recent contribution to the family of bottom-up 
approaches is the dissipative particle dynamics (DPD) 
method introduced by Hoogerbrugge and Koelman in 
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1992 ||. Although in the original formulation of DPD 
time was discrete and space continuous, a more recent re- 
interpretation asserts that this model is in fact a finite- 
difference approximation to the 'true' DPD, which is de- 
fined by a set of continuous time Langevin equations with 
momentum conservation between the dissipative parti- 
cles ||. Successful applications of the technique have 
been made to colloidal suspensions |10[], p olymer solu- 
tions (TTJ and binary immiscible fluids pl2"|. For specific 
applications where comparison is possible, this algorithm 
is orders of magnitude faster than MD p3| . The basic 
elements of the DPD scheme are particles that represent 
rather ill-defined 'mesoscopic' quantities of the underly- 
ing molecular fluid. These dissipative particles are stipu- 
lated to evolve in the same way that MD particles do, but 
with different inter-particle forces: since the DPD parti- 
cles are pictured to have internal degrees of freedom, the 
forces between them have both a fluctuating and a dissi- 
pative component in addition to the conservative forces 
that are present at the MD level. Newton's third law 
is still satisfied, however, and consequently momentum 
conservation together with mass conservation produce 
hydrodynamic behavior at the macroscopic level. 

Dissipative particle dynamics has been shown to pro- 
duce the correct macroscopic (continuum) theory; that is, 
for a one-component DPD fluid, the Navier-Stokes equa- 
tions emerge in the large scale limit, and the fluid viscos- 
ity can be computed |l4|,[l5|]. However, even though dis- 
sipative particles have generally been viewed as clusters 
of molecules, no attempt has been made to link DPD to 
the underlying microscopic dynamics, and DPD thus re- 
mains a foundationless algorithm, as is that of the hydro- 
dynamic lattice gas and a fortiori the lattice-Boltzmann 
method. It is the principal purpose of the present paper 
to provide an atomistic foundation for dissipative par- 
ticle dynamics. Among the numerous benefits gained 
by achieving this, we are then able to provide a precise 
definition of the term 'mcsoscalc', to relate the hitherto 
purely phenomcnological parameters in the algorithm to 
underlying molecular interactions, and thereby to formu- 
late DPD simulations for specific physicochemical sys- 
tems, defined in terms of their molecular constituents. 
The DPD that we derive is a representation of the under- 
lying MD. Consequently, to the extent that the approxi- 
mations made are valid, the DPD and MD will have the 
same hydrodynamic descriptions, and no separate kinetic 
theory for, say, the DPD viscosity will be needed once it 
is known for the MD system. Since the MD degrees of 
freedom will be integrated out in our approach the MD 
viscosity will appear in the DPD model as a parameter 
that may be tuned freely. 

In our approach, the 'dissipative particles' (DP) are de- 
fined in terms of appropriate weight functions that sam- 
ple portions of the underlying conservative MD parti- 
cles, and the forces between the dissipative particles are 
obtained from the hydrodynamic description of the MD 



system: the microscopic conservation laws carry over di- 
rectly to the DPD, and the hydrodynamic behavior of 
MD is thus reproduced by the DPD, albeit at a coarser 
scale. The mesoscopic (coarse-grained) scale of the DPD 
can be precisely specified in terms of the MD interac- 
tions. The size of the dissipative particles, as specified 
by the number of MD particles within them, furnishes 
the meaning of the term 'mesoscopic' in the present con- 
text. Since this size is a freely tunable parameter of the 
model, the resulting DPD introduces a general proce- 
dure for simulating microscopic systems at any conve- 
nient scale of coarse graining, provided that the forces 
between the dissipative particles are known. When a hy- 
drodynamic description of the underlying particles can 
be found, these forces follow directly; in cases where this 
is not possible, the forces between dissipative particles 
must be supplemented with the additional components 
of the physical description that enter on the mesoscopic 
level. 

The DPD model which we derive from molecular dy- 
namics is formally similar to conventional, albeit foun- 
dationless, DPD |l4| . The interactions are pairwise and 
conserve mass and momentum, as well as energy ]l6| , pT[ . 
Just as the forces conventionally used to define DPD have 
conservative, dissipative and fluctuating components, so 
too do the forces in the present case. In the present 
model, the role of the conservative force is played by 
the pressure forces. However, while conventional dissi- 
pative particles possess spherical symmetry and experi- 
ence interactions mediated by purely central forces, our 
dissipative particles are defined as space-filling cells on a 
Voronoi lattice whose forces have both central and tan- 
gential components. These features are shared with a 
model studied by Espanol |l8[ ). This model links DPD 
to smoothed particle hydrodynamics fill and defines the 
DPD forces by hydrodynamic considerations in a way 
analogous to earlier DPD models. Espanol et al. |2Cf l 
have also carried out MD simulations with a superposed 
Voronoi mesh in order to measure the coarse grained 
inter-DP forces. 

While conventional DPD defines dissipative particle 
masses to be constant, this feature is not preserved in our 
new model. In our first publication on this theory [^l[ , we 
stated that, while the dissipative particle masses fluctu- 
ate due to the motion of MD particles across their bound- 
aries, the average masses should be constant. In fact, 
the DP-masses vary due to distortions of the Voronoi 
cells, and this feature is now properly incorporated in 
the model. 

We follow two distinct routes to obtain the fluctuation- 
dissipation relations that give the magnitude of the ther- 
mal forces. The first route follows the conventional path 
which makes use of a Fokker-Planck equation ||. We 
show that the DPD system is described in an approx- 
imate sense by the isothermal-isobaric ensemble. The 
second route is based on the theory of fluctuating hy- 
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drodynamics and it is argued that this approach corre- 
sponds to the statistical mechanics of the grand canon- 
ical ensemble. Both routes lead to the same result for 
the fluctuating forces and simulations confirm that, with 
the use of these forces, the measured DP temperature is 
equal to the MD temperature which is provided as input. 
This is an important finding in the present context as the 
most significant approximations we have made underlie 
the derivation of the thermal forces. 



= ^2fk{*i)e t , (3) 

i 

where x, and v, are the position and velocity of the 
ith MD particle, which are all assumed to have identi- 
cal masses m, Pfc is the momentum of the fcth DP and 
VMD(rij) is the potential energy of the MD particle pair 
ij, separated a distance r^. The particle energy ej thus 
contains both the kinetic and a potential term. The kine- 
matic condition 



II. COARSE-GRAINING MOLECULAR 
DYNAMICS: FROM MICRO TO MESOSCALE 

The essential idea motivating our definition of meso- 
scopic dissipative particles is to specify them as clus- 
ters of MD particles in such a way that the MD par- 
ticles themselves remain unaffected while all being repre- 
sented by the dissipative particles. The independence of 
the molecular dynamics from the superimposed coarse- 
grained dissipative particle dynamics implies that the 
MD particles are able to move between the dissipative 
particles. The stipulation that all MD particles must be 
fully represented by the DP's implies that while the mass, 
momentum and energy of a single MD particle may be 
shared between DP's, the sum of the shared components 
must always equal the mass and momentum of the MD 
particle. 



A. Definitions 

Full representation of all the MD particles can be 
achieved in a general way by introducing a sampling func- 
tion 



r fc = U fe = P fc /M fe 



(4) 
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g(x - rfc) 



(1) 



where the positions and ri define the DP centers, x is 
an arbitrary position and s(x) is some localized function, 
ft will prove convenient to choose it as a Gaussian 



s(x) = exp (— x 2 /a 2 ) 



(2) 



where the distance a sets the scale of the sampling func- 
tion, although this choice is not necessary. The mass, 
momentum and internal energy E of the fcth DP are then 
defined as 



M, 



m, 



Pfc = X! /fc(xi)mvj, 

i 

X -M k XJ 2 + E k = Y, fk{*i) I \mv 2 + V MD {n 3 ) J 



completes the definition of our dissipative particle dy- 
namics. 

It is generally true that mass and momentum conserva- 
tion suffice to produce hydrodynamic behavior. However, 
the equations expressing these conservation laws contain 
the fluid pressure. In order to get the fluid pressure a 
thermodynamic description of the system is needed. This 
produces an equation of state, which closes the system 
of hydrodynamic equations. Any thermodynamic poten- 
tial may be used to obtain the equation of state. In the 
present case we shall take this potential to be the internal 
energy E k of the dissipative particles, and we shall obtain 
the equations of motion for the DP mass, momentum and 
energy. Note that the internal energy would also have to 
be computed if a free energy had been chosen for the ther- 
modynamic description. For this reason it is not possible 
to complete the hydrodynamic description without tak- 
ing the energy flow into account. As a byproduct of this 
the present DPD also contains a description of the heat 
flow and corresponds to the recently introduced DPD 
with energy conservation |Tq , p!7| . Espanol previously in- 
troduced an angular momentum variable describing the 
dynamics of extended particles fl8|| : this is needed when 
forces are non-central in order to avoid dissipation of en- 
ergy in a rigid rotation of the fluid. Angular momentum 
could be included on the same footing as momentum in 
the following developments. However for reasons both 
of space and conceptual economy we shall omit it in the 
present context, even though it is probably important in 
applications where hydrodynamic precision is important. 
In the following sections, we shall use the notation r, M, 
P and E with the indices k ,1 ,m and n to denote DP's 
while we shall use x, m, v and e with the indices i and j 
to denote MD particles. 



B. Equations of motion for the dissipative particles 
based on a microscopic description 

The fact that all the MD particles are represented at 
all instants in the coarse-grained scheme is guaranteed by 
the normalization condition J^ fc /fc(x) = 1. This implies 
directly that 
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E M * = E m 

k i 
k i 



E« 



(5) 



thus with mass, momentum and energy conserved at the 
MD level, these quantities are also conserved at the DP 
level. In order to derive the equations of motion for dissi- 
pative particle dynamics we now take the time derivatives 
of Eqs. (g). This gives 



dM fc 
dt 

dP fe 



dt 
dE tot 



= E (/fe( x O mv i + /fe(x;)F-, 

i 

(]/ = E (A(xi)e< + fk(x-i)e, 



(6) 
(7) 
(8) 



where d/di is the substantial derivative and Fj = mVj is 
the force on particle z. 

The Gaussian form of s implies that 
s(x) = — (2/a 2 )x • xs(x). This makes it possible to write 



/fc(Xj) = /fcz(xi)(v- • r k i + x- • U fei ) 



(9) 



where the overlap function is defined as /fcz(x) = 
(2/a 2 )/ fc (x)/,(x), r fci = (r fe - r,) and U fci = (U fc - U,), 
and we have rearranged terms so as to get them in terms 
of the centered variables 



v,- = v, ; - 



(Ufc + up 
2 

(r fe + r,) 



(10) 



Before we proceed with the derivation of the equations 
of motion it is instructive to work out the actual forms of 
/^(x) and /fcz(x) in the case of only two particles k and 
I. Using the Gaussian choice of s we immediately get 



Zfc(x) = 



1 



1 + [exp ((x - (r fc + r,)/2) ■ {v kl )/{^))\ 



The overlap function similarly follows: 



hi (x) = ^ cosh 



rfc + ri 



fk(r) 



A 




A 




a 2 /\r,- r.l 
k 1 

FIG. 1. The overlap region between two Voronoi cells is 
shown in grey. The sampling function /fc(r) is shown in the 
top graph and the overlap function fki{r) = (2/a 2 )fk(r)fi (r) 
in the bottom graph. The width of the overlap region is 
a 2 /|rfc — vi\ and its length is denoted by /. 

These two functions are shown in Fig.|l|. Note that the 
scale of the overlap region is not a but a 2 /|r/j — r/|. Dis- 
sipative particle interactions only take place where the 
overlap function is non-zero. This happens along the di- 
viding line which is equally far from the two particles. 
The contours of non-zero fki thus define a Voronoi lat- 
tice with lattice segments of length lu- This Voronoi 
construction is shown in Fig. in which MD particles 
in the overlap region defined by fki > 0.1, are shown, 
though presently not actually simulated as dynamic en- 
tities. The volume of the Voronoi cells will in general vary 
under the dynamics. However, even with arbitrary dis- 
sipative particle motion the cell volumes will approach 
zero only exceptionally, and even then the identities of 
the DP particles will be preserved so that they subse- 
quently re-emerge. 



1. Mass equation 



The mass equation (g) takes the form 
dM fc 



dt 



E Mm 



(11) where 



M; 



fei 



^2 /fc;(x 4 )m(v- ■ Y U + X- • Ufci) 



(13) 



(14) 



The term will be shown to be negligible within our ap- 
proximations. The x^-Ufc;-term however describes the ge- 
ometric effect that the Voronoi cells do not conserve their 
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volume: The relative motion of the DP centers causes the 
cell boundaries to change their orientation. We will re- 
turn to give this 'boundary twisting' term a quantitative 
content when the equations of motion are averaged-an ef- 
fect which was overlooked in our first publication of this 
theory pit] where it was stated that (Mm) = 0. 



f 



• t • f \ /. V 



.■<s*5«-i. 



.* 



ij ij 



Ti7 



(17) 



where Axy = x,; — Xj, we have Taylor expanded /fc(x) 
around Xj and used a result similar to Eq. (|^) to evaluate 
V/fc(x). In passing from the third to the fourth line in 
the above equations we have moved the first term on the 
right hand side to the left hand side and divided by two. 
Now, if we group the last term above with the Tki term in 
Eq. ( |l5| ) , make use of Eq. ( |To| ) , and do some rearranging 
of terms we get 

+ £/ feZ (x i )n;-r fct 

li 

+ £/ w (x i )mv^<.U fci (18) 



FIG. 2. The Voronoi lattice defined by the dissipative 
particle positions r^. The grey dots which represent the un- 
derlying MD particles are drawn only in the overlap region. 



2. Momentum equation 



The momentum equation (Q) takes the form 



£/*(x0f< 



(15) 



We can write the force as F^ = mg + Y^j F ij where the 
first term is an external force and the second term is the 
internal force caused by all the other particles. Newton's 
third law then takes the form Fjj = — Fjj. The last term 
in Eq. (Il5]) may then be rewritten as 



£ /*(x*)F< = Mkg + £ fk(xi)Fi 



(16) 



where 



where we have used the relation Mk = Yli Mki and de- 
fined the general momentum-flux tensor 



n, 



i^FyAxi 



(19) 



This tensor is the momentum analogue of the mass-flux 
vector mvj. The prime indicates that the velocities on 
the right hand side are those defined in Eq. ( |l0| ) . The ten- 
sor II,; describes both the momentum that the particle 
carries around through its own motion and the momen- 
tum exchanged by inter-particle forces. It may be arrived 
at by considering the momentum transport across imag- 
inary cross sections of the volume in which the particle 
is located. 



3. Energy equation 

In order to get the microscopic energy equation of mo- 
tion we proceed as with the mass and momentum equa- 
tions and the two terms that appear on the right hand 
side of Eq. (§). 

Taking Vmd to be a central potential and using 
the relations VVMD(fy) = ^Mi)( r y ) e " = ~ F " an d 
VMD(nj) = V^ D (rij)eij -Vy = -Fij ■ 



where v, 



Vi — Vj we get the time rate of change of the particle 
energy 
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(20) 



This gives the first term of Eq. (^) in the form 
E fk(*i)i = Pk ■ g + \ E / fe (x i )F ij • (v 4 + VJ ) . (21) 

i i^tj 

The last term of this equation is odd under the exchange 
i <-» j and exactly the same manipulations as in Eq. Jl7| ) 
may be used to give 

^/ fc (x. ( )e = P fe ■ g 

i 

+ E /w( X i)T F ij ' ( V ^ + V i) Ax ij ' V kl 

= P fc -g+^/ H (xi)(jFy-(v' i + v<) 

+ ~Fy • Ufc + U ^ Axy • r fci (22) 

where for later purposes we have used Eqs. ( [Toj ) to get the 
last equation. The last term of Eq. (|J) is easily written 
down using Eq. (|^). This gives 

E A-( x iK = E /fe'( x 0( v i ' r fc; + x i ' UfciK . (23) 

i ii 

As previously we write the particle velocities in terms of 
v.; . The corresponding expression for the particle energy 
is e 4 = e[ + rnvj • (U fc + U,)/2 + (l/2)m((U fe + U,)/2) 2 
where the prime in e[ denotes that the particle velocity 
is v£ rather than Vj. Equation (p3|) may then be written 



+ E /fci( x [ e i v i + • — 2 I ' r -< 

^/ H (x 4 )^-U fei . (24) 



Combining this equation with Eq. (|22j) we obtain 



E/HW(^+^'('^||x; U t( . ,25) 



where the momentum-flux tensor is defined in Eq. ( |l9| ) 
and we have identified the energy-flux vector associated 
with a particle i 



z£*< 



(v< 



Vj)AXy 



(26) 



Again the prime denotes that the velocities are V { rather 
than Vj . To get the internal energy E^ instead of E^ we 
note that d(P 2 k /2M k )/dt = U fe • P fe - (l/2)M fc U|. Using 
this relation, the momentum equation Eq. ( |l8l ) , as well as 
the substitution (U fc + U/)/2 = U fc - Ujm/2 in Eq. (§|), 
followed by some rearrangement of the M^; terms we find 
that 



dt \2 



Ui 



u, 



E/«( 



E /«(*)( J* -n< .> 

, Uki 



X, ) I f ; - ///V; ■ — - I • U fe ; 



This equation has a natural physical interpretation. 
The first term represents the translational kinetic energy 
of the DP as a whole. The remaining terms represent 
the internal energy E k . This is a purely thermodynamic 
quantity which cannot depend on the overall velocity of 
the DP, i.e. it must be Galilean invariant. This is eas- 
ily checked as the relevant terms all depend on velocity 
differences only. 

The Mki term represents the kinetic energy received 
through mass exchange with neighboring DPs. As will 
become evident when we turn to the averaged descrip- 
tion, the term involving the momentum and energy fluxes 
represents the work done on the DP by its neighbors and 
the heat conducted from them. The e^-term represents 
the energy received by the DP due to the same 'bound- 
ary twisting' effect that was found in the mass equation. 
Upon averaging, the last term proportional to will be 
shown to be relatively small since (vJ) = in our approx- 
imations. This is true also in the mass and momentum 
equations. Equations (fTi]), (18) and ( |27| ) have the coarse 
grained form that will remain in the final DPD equations. 
Note, however, that they retain the full microscopic in- 
formation about the MD system, and for that reason they 
are time-reversible. Equation ( p"8| ) for instance contains 
only terms of even order in the velocity. In the next sec- 
tion terms of odd order will appear when this equation 
is averaged. 

It can be seen that the rate of change of momentum 
in Eq. (^) is given as a sum of separate pairwise con- 
tributions from the other particles, and that these terms 
are all odd under the exchange I <-> k. Thus the parti- 
cles interact in a pairwise fashion and individually fulfill 
Newton's third law; in other words, momentum conser- 
vation is again explicitly upheld. The same symmetries 
hold for the mass conservation equation ( p^ ) and energy 
equation (Eq). 
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III. DERIVATION OF DISSIPATIVE PARTICLE 
DYNAMICS: AVERAGE AND FLUCTUATING 
FORCES 

We can now investigate the average and fluctuating 
parts of Eqs. (§?]), (|l|) and Q. In so doing we shall 
need to draw on a hydrodynamic description of the un- 
derlying molecular dynamics and construct a statistical 
mechanical description of our dissipative particle dynam- 
ics. For concreteness we shall take the hydrodynamic de- 
scription of the MD system in question to be that of a 
simple Newtonian fluid [Q. This is known to be a good 
description for MD fluids based on Lennard- Jones or hard 
sphere potentials, particularly in three dimensions ||. 
Here we shall carry out the analysis for systems in two 
spatial dimensions; the generalization to three dimen- 
sions is straight forward, the main difference being of 
a practical nature as the Voronoi construction becomes 
more involved. 

We shall begin by specifying a scale separation between 
the dissipative particles and the molecular dynamics par- 
ticles by assuming that 



Xi - x,- << \r k - n\ 



(28) 



where x^ and Xj denote the positions of neighbouring MD 
particles. Such a scale separation is in general necessary 
in order for the coarse-graining procedure to be physi- 
cally meaningful. Although for the most part in this pa- 
per we are thinking of the molecular interactions as being 
mediated by short-range forces such as those of Lennard- 
Jones type, a local description of the interactions will still 
be valid for the case of long-range Coulomb interactions 
in an electrostatically neutral system, provided that the 
screening length is shorter than the width of the over- 
lap region between the dissipative particles. Indeed, as 
we shall show here, the result of doing a local averag- 
ing is that the original Newtonian equations of motion 
for the MD system become a set of Langevin equations 
for the dissipative particle dynamics. These Langevin 
equations admit an associated Fokker-Planck equation. 
An associated fluctuation-dissipation relation relates the 
amplitude of the Langevin force to the temperature and 
damping in the system. 



A. Definition of ensemble averages 

With the mesoscopic variables now available, we need 
to define the correct average corresponding to a dynam- 
ical state of the system. Many MD configurations are 
consistent with a given value of the set {r^, Mk, Ufe, Ek}, 
and averages are computed by means of an ensemble of 
systems with common instantaneous values of the set 
{rfc, Mfc, Ufe, Ek}- This means that only the time deriva- 
tives of the set {r^, Mk, Ek}, i.e. the forces, have a 



fluctuating part. In the end of our development approxi- 
mate distributions for UVs and E^s will follow from the 
derived Fokker-Planck equations. These distributions re- 
fer to the larger equilibrium ensemble that contains all 
fluctuations in {r^, Mk, Ufc, Ek}- 

It is necessary, to compute the average MD particle 
velocity (v) between dissipative particle centers, given 
{rfc, Mk, Ufc, Ek}- This velocity depends on all neighbor- 
ing dissipative particle velocities. However, for simplicity 
we shall only employ a "nearest neighbor" approxima- 
tion, which consists in assuming that (v) interpolates lin- 
early between the two nearest dissipative particles. This 
approximation is of the same nature as the approxima- 
tion used in the Newtonian fluid stress-strain relation 
which is linear in the velocity gradient. This implies that 
in the overlap region between dissipative particles fc and 
I 



<v') = (v')(x) = 



Ykl 



U 



(29) 



kl 



where the primes are defined in Eqs. (10) and Tki = 
|r fc -rj|. 

A preliminary mathematical observation is useful in 
splitting the equations of motion into average and fluc- 
tuating parts. Let r(x) be an arbitrary, slowly varying 
function on the a 2 jrki scale. Then we shall employ the 
approximation corresponding to a linear interpolation be- 
tween DP centers, that r(x) = (l/2)(rk + rj) where x is 
a position in the overlap region between DP k and 1 and 
Tk and Ti are values of the function r associated with the 
DP centers k and 1 respectively. 




kl 



1 



J kl 



k 




FIG. 3. Two interacting Voronoi cells. The length of the 
intersection between DP's k and 1 is Iki, the shift from the 
center of the intersection between r^i and Iki is Lki {Lki = 
when Tki intersects lki in the middle) and the unit vector iki 
is normal to e^. The coordinate system x-y used for the 
integration has its origin on the intersection. 
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Then 



Ei/w(xi)r(x) fa J dx dy Pk ^ Pl fki{x) Tk ^ n 



hi Pk + Pir k + n 
2a 2 2 2 
^ Pk + Pir k + n 
r k i 2 2 ' 



da;' cosh (x'rki/a ) 



(30) 



where Pfc + Pi is the MD particle number density and we 
have used the identity tanh'(x) = cosh - (x). We will 
also need the first moment in x' 

^/fci(xi)x-r(xi) pa y dxdy^^^-f M (x)x'^-^- 

i 

J dx dycosh -2 (^r) J/ifci 
(31) 



1 Pfc + Pi r fc + r, 



2a 2 2 I 
iki r Pk+ Pi r k + n . 

— J^kl lH 

2r u 2 2 

where the unit vectors e k i = r k i/ru and i k i are shown 
in Fig. ||, we have used the fact that the integral 
over xe k i cosh -2 ... vanishes since the integrand is odd, 
and the last equation follows by the substitution x — > 
(a 2 /r k i)x. In contrast to the vector e k i the vector i k i is 
even under the exchange k <-> i, as is Lfcz. This is a mat- 
ter of definition only as it would be equally permissible 
to let iki and Lj-i be odd under this exchange. However, 
it is important for the symmetry properties of the fluxes 
that hi and L k i have the same symmetry under k <-> I. 



B. The mass conservation equation 

Taking the average of Eq. (|l4|), we observe that the 
first term vanishes if Eq. (^9|) is used, and the second 
term follows directly from Eq. (|3l|). We thus obtain 



M fe = ^((M H ) + id kl ) 



where 



(M, 



kl I 



E 



/ fc im(xi)(x-) • XJ k i 



hi T Pk + 

J-ikl 

2r kt 2 



Pi 



(32) 



hi ■ U A 



(33) 



and Mti = Mki — (Mm)- The finite value of (Mm) is 
caused by the relative DP motion perpendicular to e^;. 
This is a geometric effect intrinsic to the Voronoi lat- 
tice. When particles move the Voronoi boundaries change 
their orientation, and this boundary twisting causes mass 
to be transferred between DP's. This mass variation will 
be visible in the energy flux, though not in the momen- 
tum flux. It will later be shown that the effect of mass 
fluctuations in the momentum and energy equations may 
be absorbed in the force and heat flux fluctuations. 



C. The momentum conservation equation 



Using Eq. (^) we may split Eq. (|l^) into average and 
fluctuating parts to get 



dP fc 
dt 



M fc g 

i u 
+ Y, /*' ( x *Mv; ; x' t ) • U kl + X) F W , (34) 



where the fluctuating force or, equivalently, the momen- 
tum flux is 

FM=X;/M(xi)[(IIi-(IIi))-r H 



+ ™(v-x- - (vjxl)) • Ufci] 



(35) 



Note that by definition Fjfc = —F k i- The fact that we 
have absorbed mass fluctuations with the fluctuations in 
Ffcj deserves a comment. In general force fluctuations 
will cause mass fluctuations, which in turn will couple 
back to cause momentum fluctuations. The time scale 
over which this will happen is t v = r^/r], where r\ is the 
dynamic viscosity of the MD system. This is the time 
it takes for a velocity perturbation to decay over a dis- 
tance of r k i- Perturbations mediated by the pressure, i.e. 
sound waves, will have a shorter time. In the sequel we 
shall need to make the assumption that the forces are 
Markovian, and it is clear that this assumption may only 
be valid on time scales larger than t v . Since the time 
scale of a hydrodynamic perturbation of size I, say, is 
also given as l 2 /n this restriction implies the scale sepa- 
ration requirement << I 2 , consistent with the scale 
r k i being mesoscopic. 

Since (Ilj) is in general dissipative in nature, Eq. ( |3^ ) 
and its mass- and energy analogue will be referred to as 
DPD1. It is at the point of taking the average in Eq. ( |3^ ) 
that time reversibility is lost. Note, however, that we do 
not claim to treat the introduction of irreversibility into 
the problem in a mathematically rigorous way. This is 
a very difficult problem in general which so far has only 
been realized by rigorous methods in the case of some 
very simple dynamical systems with well defined ergodic 
properties p2]-p4|. We shall instead use the constitu- 
tive relation for a Newtonian fluid which, as noted ear- 
lier, is an emergent property of Lennard-Jones and hard 
sphere MD systems, to give Eq. (34) a concrete content. 
The momentum-flux tensor then has the following simple 
form 
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p(Hi) — mpvv + Ip — ?7(Vv + (Vv) 



(36) 



where p is the pressure and v the average velocity of the 
MD fluid, T denotes the transpose and I is the identity 
tensor In the above equation we have for simplicity 
assumed that the bulk viscosity ( = (2/d)rj where d is 
the space dimension 2. The modifications to include an 
independent £ are completely straight forward. 

Using the assumption of linear interpolation (Eq. (p9|)), 
the advective term pvv vanishes in the frame of reference 
of the overlap region since there v'«0. The velocity gra- 
dients in Eq. ([56j) may be evaluated using Eq. (p9[); the 
result is 

Vv + (Vv) T - — {e kl XJ kl + U w e«) . (37) 

Note further that J2i hi is in fact a surface integral 
over the DP surface. Consequently 



( 



hl^kl9k — 



(38) 



for any function g k that does not depend on I. In par- 
ticular we have Y^ihl^klipk +j>j)/2 = -Y^kiekiPki/Z 
where p k i — pk—pi- Combining Eqs. j3fj), (|30|) and (37) 
Eq. (U) then takes the form 

Mfcg + 2^{M U ) = 



df 

i 



i 

/./ 1 ^7T e ki + — (Ufei + (Ufc( • e k i)e k i 
2 r fc; 



fei 



(39) 



where we have assumed that the pressure p, as well as the 
average velocity, interpolates linearly between DP cen- 
ters, and we have omitted the (v£x£) « term. Note that 
all terms except the gravity term on the right hand side of 
Eq. ([39|) are odd when k «-> I. This shows that Newton's 
third law is unaffected by the approximations made and 
that momentum conservation holds exactly. The same 
statements can be made for the mass equation and the 
energy equation. The pressure will eventually follow from 
an equation of state of the form p k = p(E k , V k , M k ) where 
Vk is the volume and M k is the mass of DP k. 



D. The energy conservation equation 



Splitting Eq. (27) into an average and a fluctuating 
part gives 



u 



kl 



r k i 



I v 

* kl ■ r q k i ■ 



(40) 



where we have defined 

Qkl = fkl(*i){ J 'et ~ (J«>) ' r kl + ( 



U 



kl in tt 



(41) 



i.e. the fluctuations in the heat flux also contains the 
energy fluctuations caused by mass fluctuations. This is 
like the momentum case. 

Note that in taking the average in Eq. ( ftO| ) the II • 
product presents no problem as XJ k i is kept fixed under 
this average. If we had averaged over different values of 
U k i the product of velocities in Il-Ujy would have caused 
difficulties. Equation ( ETj ) is the third component in the 
description at the DPD1 level. 

The average of the energy flux vector J e is taken to 
have the general form M 



AVT 



(42) 



where a = II — pvv is the stress tensor, and A the ther- 
mal conductivity and T the local temperature. Note 
that in Eq. (|27]) only 3' t appears. Since v' w we have 
(J' e ) = AVT. Averaging of Eq. © gives 



Tki 



E, jPk+Pl V /t t . f-r t \ A Ufci 

Ilk — - — e u ( U ki + ( U kl ■ efcj )eki ) ■ ——- 
v 2 r u I 2 



El/w x /U/d\ hi . T / -E& -E; 
^ 2 V 2 / 4r fei V Vfc Vi 

f fcZ ■ -r — h ftZ • 



(43) 



where Tjy = T& — T; is the temperature difference be- 
tween DP's k and /, and we have used linear interpola- 
tion to write (ei) = (\/2)(E k /V k + Ei/V{). The first term 
above describes the heat flux according to Fourier's law. 
The next non-fluctuating terms, which are multiplied by 
Ufcz/2 represent the (rate of) work done by the interpar- 
ticle forces, and the F k [ term represents the work done 
by the fluctuating force. 

As has been pointed out by Avalos et al and Espanol 
p6| , p7t the work done by Fki has the effect that it in- 
creases the thermal motion of the DP's at the expense 
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of a reduction in Ek . This is the case here as well since 
the above Fki • ^Ju term always has a positive average 
due to the positive correlation between the force and the 
velocity increments. 

Equation ( pL3| ) is identical in form to the energy equa- 
tion postulated by Avalos and Mackie ITS], save for the 
fact that here the conservative force ((pfc + pi)/2)eki ■ 
Ufei/2 (which sums to zero under J2k) * s present. The 
pressure forces in the present case correspond to the con- 
servative forces in conventional DPD-it will be observed 
that they are both derived from a potential. However, 
while the conservative force in conventional DPD must 
be thought to be carried by some field external to the 
particles, the pressure force in our model has its origin 
within the particles themselves. There is also a small 
difference between the present form of Fourier's law and 
the description of thermal conduction employed by Ava- 
los and Mackie. While the heat flux here is taken to be 
linear in differences in T, Avalos and Mackie use a flux 
linear in differences in (1/T). As both transport laws are 
approximations valid to lowest order in differences in T, 
they should be considered equivalent. 

With the internal energy variable at hand it is possible 
to update the pressure and temperature T of the DP's 
provided an equation of state for the underlying MD sys- 
tem is assumed, and written in the form P — P{E, V, m) 
and T — T(E, V, m). For an ideal gas these are the well 
known relations PV = (2/d)E and k B T = (2/d)mE. 

Note that we only need the average evolution of the 
pressure and temperature. The fluctuations of p are al- 
ready contained in Fki and the effect of temperature fluc- 
tuations is contained within q\i ■ 

At this point we may compare the forces arising in 
the present model to those used in conventional DPD. 
In conventional DPD the forces are pairwise and act in 
a direction parallel to e^i , with a conservative part that 
depends only on rjy and a dissipative part proportional 
to (Ufc; • ejy)e&; m^lH^]. The forces in our new version 
of DPD are pairwise too. The analog of the conserva- 
tive force, lkl(pkl/2)&kU is central and its r dependence 
is given by the Voronoi lattice. When there is no overlap 
Iki between dissipative particles their forces vanish. (A 
cut-off distance, beyond which no physical interactions 
are permitted, was also present in the earlier versions of 
DPD-see, for example, Ref. (|]-where it was introduced 
to simplify the numerical treatment.) Due to the exis- 
tence of an overlap region in our model, the dissipative 
force has both a component parallel to and a com- 
ponent parallel to the relative velocity XJki- However, 
due to the linear nature of the stress-strain relation in 
the Newtonian MD fluid studied here, this force has the 
same simple linear velocity dependence that has been 
postulated in the literature. 

The friction coefficient is simply the viscosity r\ of the 
underlying fluid times the geometric ratio Iki jtki ■ As has 
been pointed out both in the context of DPD |Q and 



elsewhere, the viscosity is generally not proportional to a 
friction coefficient between the particles. After all, con- 
servative systems like MD are also described by a viscos- 
ity. Generally the viscosity will be caused by the com- 
bined effect of particle interaction (dissipation, if any) 
and the momentum transfer caused by particle motion. 
The latter contribution is proportional to the mean free 
path. The fact that the MD viscosity r], the DPD vis- 
cosity and the friction coefficient are one and the same 
therefore implies that the mean free path effectively van- 
ishes. This is consistent with the space filling nature of 



the particles. See Sec. VI B for a further discussion of the 
zero viscosity limit. 

Note that constitutive relations like Eqs. (E36|) and (|4S 



are usually regarded as components of a top-down or 
macroscopic description of a fluid. However, any bottom- 
up mesoscopic description necessarily relies on the use of 
some kind of averaging procedure; in the present con- 
text, these relations represent a natural and convenient 
although by no means a necessary choice of average. The 
derivation of emergent constitutive relations is itself part 
of the programme of non-equilibrium statistical mechan- 
ics (kinetic theory), which provides a link between the 
microscopic and the macroscopic levels. However, as 
noted above, no general and rigorous procedure for de- 
riving such relations has hitherto been realised; in the 
present theoretical treatment, such assumed constitutive 
relations are therefore a necessary input in the linking of 
the microscopic and mesoscopic levels. 



IV. STATISTICAL MECHANICS OF 
DISSIPATIVE PARTICLE DYNAMICS 

In this section we discuss the statistical properties of 
the DP's with the particular aim of obtaining the magni- 
tudes of Fki and q\i . We shall follow two distinct routes 
that lead to the same result for these quantities, one 
based on the conventional Fokker-Planck description of 
DPD ||l6|, and one based on Landau's and Lifshitz's fluc- 
tuating hydrodynamics |jj . 

It is not straightforward to obtain a general statistical 
mechanical description of the DP-system. The reason is 
that the DP's, which exchange mass, momentum, energy 
and volume, are not captured by any standard statistical 
ensemble. For the grand canonical ensemble, the sys- 
tem in question is defined as the matter within a fixed 
volume, and in the case of a the isobaric ensemble the 
particle number is fixed. Neither of these requirements 
hold for a DP in general. 

A system which exchanges mass, momentum, energy 
and volume without any further restrictions will gener- 
ally be ill-defined as it will lose its identity in the course 
of time. The DP's of course remain well-defined by virtue 
of the coupling between the momentum and volume vari- 
ables: The DP volumes are defined by the positions of the 
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DP-centers and the DP-momenta govern the motion of 
the DP-centers. Hence the quantities that are exchanged 
with the surroundings are not independent and the en- 
semble must be constructed accordingly. 

However, for present purposes we shall leave aside the 
interesting challenge of designing the statistical mechan- 
ical properties of such an ensemble, and derive the mag- 
nitude of Fki and q k i from two different approximations. 
The approximations are both justifiable from the as- 
sumption that Ffci and qu have a negligible correlation 
time. It follows that their properties may be obtained 
from the DP behavior on such short time scales that the 
DP-centers may be assumed fixed in space. As a result, 
we may take either the DP volume or the system of MD- 
particles fixed for the relevant duration of time. Hence 
for the purpose of getting Fki and qu we may use cither 
the isobaric ensemble, applied to the DP system, or the 
grand canonical ensemble, applied to the MD system. 
We shall find the same results from either route. The 
analysis of the DP system using the isobaric ensemble 
follows the standard procedure using the Fokker-Planck 
equation, and the result for the equilibrium distribution 
is only valid in the short time limit. The analysis of the 
MD system corresponding to the grand canonical ensem- 
ble could be conducted along the similar lines. However, 
it is also possible to obtain the magnitude of Ffc; and 
qki directly from the theory of fluctuating hydrodynam- 
ics since this theory is derived from coarse-graining the 
fluid onto a grid. The pertinent fluid velocity and stress 
fields thus result from averages over fixed volumes associ- 
ated with the grid points: Since mass flows freely between 
these volumes the appropriate ensemble is thus the grand 
canonical one. 



and our analysis parallels that of Avalos and Mackie [[16| . 
However, the fact that the conservative part of the con- 
ventional DP forces is here replaced by the pressure and 
that the present DP's have a variable volume makes a 
separate treatment enlightening. 

The probability p(Vk, Pfc, E k ) of finding DPfc with a 
volume Vfc, momentum Pfc and internal energy E k is then 
proportional to exp^Sr/hs) where St is the entropy of 
all DP's given that the values (Vfc, Pfc, Ek) are known for 
DPfc p6|. If S' denotes the entropy of the bath we can 
write St as 



S T = S'(V T - V fc , P T - Pfc, E T - — ^ - E h ) + S k 

ZM k 



8<3' ( P 2 



dS' 

ap 



Pfc + s k 



(44) 



where the derivatives are evaluated at (Vt,Pt, Et) and 
thus characterize the bath only. Assuming that P^ van- 
ishes there is nothing in the system to give the vector 
dS' /dP a direction, and it must therefore vanish as well 
p7[ . The other derivatives give the pressure po and tem- 
perature To of the bath and we obtain 



St — S (Vtj Ptj Et) 



1 



pi 

Gfc+ Fk 



2M fe 



(45) 



where the Gibbs free energy has the standard form G k — 
E k + poV k — XbSfc. Since there is nothing special about 
DPfc it immediately follows that the the full equilibrium 
distribution has the form 



A. The isobaric ensemble 



= Z 1 (T ,p )exp I -A) ^ 



P 2 

2M i: 



Gfc 



(46) 



We consider the system of N k ^> 1 MD particles inside 
a given DPfc at a given time, say all the MD particles 
with positions that satisfy /fe(xi) > 1/2 at time to- At 
later times it will be possible to associate a certain vol- 
ume per particle with these particles, and by definition 
the system they form will exchange volume and energy 
but not mass. We consider all the remaining DP's as a 
thermodynamic bath with which DPfc is in equilibrium. 
The system defined in this way will be described by the 
Gibbs free energy and the isobaric ensemble. Due to 
the diffusive spreading of MD-particles, this system will 
only initially coincide with the DP; during this transient 
time interval, however, we may treat the DP's as sys- 
tems of fixed mass and describe them by the approxima- 
tion (M k i) = 0. The magnitudes of q and F follow in the 
form of fluctuation-dissipation relations from the Fokker- 
Planck equivalent of our Langevin equations. The math- 
ematics involved in obtaining fluctuation-dissipation re- 
lations is essentially well-known from the literature M, 



where (3o = l/(fcsTo). The temperature T k = 
(dSk/dEk)^ 1 and pressure pfc = T k (dS k /dV k ) will fluctu- 
ate around the equilibrium values To and po- The above 
distribution is analyzed by Landau and Lifshitz p7| who 
show that the fluctuations have the magnitude 



(APfc 2 ) 



VfcKs' 



(AT fe 2 ) = 



Vc v 



(47) 



where the isentropic compressibility k$ = 
— {l/V){dV/dP)s and the specific heat capacity c„ are 
both intensive quantities. Comparing our expression 
with the distribution postulated by Avalos and Mackie, 
we have replaced the Helmholtz by the Gibbs free en- 
ergy in Eq. (|46|). This is due to the fact that our DP's 
exchange volume as well as energy. 
We write the fluctuating force as 



' ki 



(48) 
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where, for reasons soon to become apparent, we have 
chosen to decompose F k i into components parallel and 
perpendicular to e^. The W's are defined as Gaussian 
random variables with the correlation function 

(W k i a (t)W nmf }(t')) = S a0 S(t - t')(6 ) (49) 

where a and f3 denote either _L or ||. The product of 5 
factors ensures that only equal vectorial components of 
the forces between a pair of DP's are correlated, while 
Newton's third law guarantees that u} k i = —i^ik- Like- 
wise the fluctuating heat flux takes the form 



q kt = A kl W kl (50) 
) without the 5 a p factor and 



where W k i satisfies Eq. ( 
energy conservation implies A k i — —A;/.. 

The force correlation function then takes the form 



(F kn (t)F lm {t')) 



nm H~ &kr< 
n ii 



S ln )S(t-t') 

+ 5kmSln)S(t 



(51) 



where we have introduced the second order tensor u) kn i m . 

It is a standard result in non-equilibrium statistical 
mechanics that a Langevin description of a dynamical 
variable y 

y = a(y) + G (52) 

where G is a delta-correlated force has an equiva- 
lent probabilistic representation in terms of the Fokker- 
Planck equation 

= -V ■ (a(y)p(y)) + IvV: (A(y)p(y)) (53) 

where V denotes derivatives with respect to y and p(y, t) 
is the probability distribution for the variable y at time 
t, (G(y, <)G(y, t')) — A8(t — t') and A is a symmetric 
tensor of rank two J28|. 

In the preceding paragraph, G denotes all the fluc- 
tuating terms in Eqs. d39| ) and ([l3|). Using the above 
definitions and (M k i) = it is a standard matter j| to 
obtain the Fokker-Planck equation 

dp 
dt 



(L + L D i S + Ldif), P 



(54) 



where 



d p k 
e k i • yj k i — 







Pi 



efer 



Pki 



L 



DIS 



J DIF 



8E k 

£'« 

kytl 

-y 

2 ^ 

k^l 

2 I 

kl 
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\dP k 



* kl - 


d 
8E k 




' t kl~ 


A— — 

ru 


d 
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Uki 




dE k 


2 



A 



\dE k 



d 

m 



(55) 



F u = (v/ r ki)(Vki + (U fci • e u )e k i), and the sum J2 k =ii 
runs over both k and I. The operator C k i is defined as in 
Ref. M: 



£ki — 



( d 







\dP k dVi 



Uh (_d_ _3_ 

2 \8E k dE l 



(56) 



The steady-state solution of Eq. ( |54| ) is already given 
by Eq. (46); following conventional procedures we can 



obtain the fluctuation-dissipation relations for lj and A 
by inserting p cq in Eq. (|5~4|). 

Apart from the tensorial nature of Wkiki the operators 
Ldis and £dif are essentially identical to those published 
earlier in conventional DPD ||l6| , [l7|] . However, the 'Liou- 
ville' operator Lq plays a somewhat different role as it 
contains the d/dE k term, corresponding to the fact that 
the pressure forces do work on the DP's to change their 
internal energy. 

While Lop cq conventionally vanishes exactly by con- 
struction of the inter-DP forces, here it vanishes only to 
order 1/N k . In order to evaluate Lop eq we need the fol- 
lowing relationship 



d 
dr k 



klGkl 



kjtl 



d 

dVi 



d 
dV k 



(57) 



which is derived by direct geometrical consideration of 
the Voronoi construction. By repeated use of Eq. (|3|) it 
is then a straightforward algebraic task to obtain 



LoP° 



4 ^ 

k^l 



hi e ki ■ U* 



dpi 
dEi 



PkiT t 



k B T k Ti 



(58) 



which docs not vanish identically. However, note that 
if we estimate Ei m NiksT we obtain dpi/dEi w 
(l/Nkfipi/ksT). Similarly we may estimate p k i and T k i 
from Eq. (E7f to obtain 



p kl T kl VAP 2 AT 2 



ksT k Ti 



ksT k Ti 



Nk/Vk 



(59) 



The last square root is an intensive quantity of the or- 
der Pq/(&bTq), as may be easily demonstrated for the 
case of an ideal gas. Since each separate quantity that 
is contained in the differences in the square brackets of 
Eq. flSq ) is of the order Po/Tq we have shown that they 
cancel up to relative order 1/N k « 1. In fact, it is not 
surprising that Langevin equations which approximate 
local gradients to first order only in the corresponding 
differences, like T k i, give rise to a Fokker-Planck descrip- 
tion that contains higher order correction terms. 
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Having shown that Lop cq vanishes to a good ap- 
proximation we may proceed to obtain the fluctuation- 
dissipation relations from the equation (Ldis + 
Ldif)p cc| = 0. It may be noted from Eq. ( |55| ) that this 
equation is satisfied if 



1 



{l kl F k J l + -u Mkl C kl )p c(i = Q 



kl 



2 f_d_ 

\dE k 



d 
dEi 



. 



Using the identity 



e k ie k i + hihi — I 

where i k i a vector normal to e k i, 
Eq. ( |60| ) implies that 



(60) 



(61) 



we may show that 



u\\ 



A 2 



2u>li ± = Ar]k B Q k i — 

2k B T k T l \— , 
r k i 



(62) 



where 



kl 



(1/2XT,- 1 + Tf 1 ). 



B. F from fluctuating hydrodynamics 

Having derived the fluctuation-dissipation relations 
from the approximation of the isobaric ensemble we now 
derive the same result from fluctuating hydrodynamics, 
which corresponds to the grand canonical ensemble. We 
shall only derive the magnitude of F k i since q follows on 
the basis of the same reasoning. 

Fluctuating hydrodynamics jjj is based on the conser- 
vation equations for mass, momentum and energy with 
the modification that the momentum and energy fluxes 
contain an additional fluctuating term. Specifically, the 
momentum flux tensor takes the form —VP + pvv + a' , 
where P is the pressure, v is the velocity field and the 
viscous stress tensor is given as 



a' = ri ( Vv + Vv 3 



• CV- 



(63) 



where s is the fluctuating component of the momentum 
flux. From the same approximations as we used in de- 
riving Eq. ([62]), i.e. a negligible correlation time for the 
fluctuating forces, Landau and Lifshitz derive 



(s(x, t) • ns(x', 0) • n) = 2k B T ( r?(l + nn) + {( - -r))nn 



6(t) 



3^ ifx,x' £ AK 
otherwise 



(64) 



where n is an arbitrary unit vector and n labels the vol- 
ume element AV^. By following the derivations presented 



by Landau and Lifshitz, it may be noted that nowhere is 
it assumed that the AT/^'s are cubic or stationary. 




FIG. 4. Pairwise Voronoi-cell interactions. Dark gray: 
The volume Vki associated with the interaction between a 
single DP-pair. The light gray region shows the volume of 
the neighboring interaction. 

By making the identifications C = {2/d)rj, n — > 
F k i = l k is ■ e kl , AV n -> Vm, (shown in Fig. |), and 
T = Q k i we may immediately write down 



(F H (*)F„ m (0)) =2 



e k ie nm )5(t) 



( S kn Si m + S km Sl n ) 



(65) 



where again the last sum of (5-factors ensures that kl and 
nm denote the same DP pair. Observing from Fig. |J that 
V k i = lufkh it now follows directly from Eq. ( |65| ) that 



(F/d (t) ■ e (0)- 



kl ~^nm (0) • i„ m ) 



= Ak B e kl l — v s(t) 

r k i 

( 5 kn 8l rn + S km Sl n ) 



(66) 



which is nothing but the momentum part of Eq. (|62j). 
That the fluctuating heat flux q produces the form of 
fluctuation-dissipation relations given in Eq. ( |62|) follows 
from a similar analysis. Thus the approximation of fixed 
DP volume V k produces the same result as the approxi- 
mation of fixed number of MD particles N k . This is due 
to the fact that both approximations are based on the 
assumption that the DP's are only considered within a 
time interval which is longer than the correlation time of 
the fluctuations but shorter than the time needed for the 
DP's to move significantly. 

The result given in Eq. ( pfj| ) was derived from the some- 
what arbitrary choice of discretizaton volume V k i ; this is 
the volume which corresponds to the segment l k i over 
which all forces have been taken as constant. It is thus 
the smallest discretization volume we may consistently 
choose. It is reassuring that Eq. (p6h also follows from 
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different choices of AV n . For example, one may readily 
check that Eq. ( |66|) is obtained if we split V k i in two along 
r k i and consider to be the sum of two independent 
forces acting on the two parts of Zjy. 

We are now in a position to quantify the av- 
erage component (E k ) = X) 



U w /2) of 

the fluctuations in the internal energy given in 
Eq. (|43l). Writing the velocity in response to F k i as 

Ufc = J2i^ k r-oo dt '^l(t')/M k , we get that {k k ) = 
E/!oo dt'(F kl (t')F kl (t)) which by Eqs. © and @ be- 
comes (E k ) = (l/M k )J23lkivkB<dki/r k i- This result is 
the same as one would have obtained applying the rules of 
Ito calculus to JJ\/ (2M k ). It yields the modified, though 
equivalent, energy equation 



E k — — ^ h k \ 



T k i 
r k i 



Pk + pi 



v_ 

rki 



II, 1 ^ efe( _ ~(^ kl + ' e kl) e kl, 



V* F' fe; • -— - 3—r]k B &kl + m 
^ 2 r ki 



Ufci 
2 

(67) 



where we have written F' k i with a prime to denote that 
it is uncorrelated with XJ k i. In a numerical implemen- 
tation this implies that F' k i must be generated from a 
different random variable than F k i, which was used to 
update Ufcz. 

The fluctuation-dissipation relations Eqs. ( |62| ) com- 
plete our theoretical description of dissipative particle 
dynamics, which has been derived by a coarse-graining 
of molecular dynamics. All the parameters and proper- 
ties of this new version of DPD are related directly to 
the underlying molecular dynamics, and properties such 
as the viscosity which are emergent from it. 




FIG. 5. The DPD temperature (energy units) averaged 
over 5000 dissipative particles as a function of time (iteration 
number in the integration scheme) , showing good convergence 
to the underlying molecular dynamics temperature which was 
set at one. This simulation provides strong support for the 
approximations used to derive the fluctuation-dissipation re- 
lations in our DPD model from molecular dynamics. 

Figure [5] shows the relaxation process towards equi- 
librium of an initially motionless system. The DP tem- 
perature is measured as (P^/ '(2M k )) for a system of DPs 
with internal energy equal to unity. The simulations were 
run for 4000 iterations of 5000 dissipative particles and 
a timestep dt — 0.0005 using an initial molecular density 
p = 5 for each DP. The molecular mass was taken to 
be m = 1, the viscosity was set at 77 = 1, the expected 
mea n free path is 0.79, and the Reynolds number (See 
Sec. |VI B| ) is Re=2.23. It is seen that the convergence of 
the DP system towards the MD temperature is good, a 
result that provides strong support for the fluctuation- 
dissipation relations of Eq. (64) . 



VI. POSSIBLE APPLICATIONS 



A. Multiscale phenomena 



V. SIMULATIONS 



While the present paper primarily deals with theoreti- 
cal developments we have carried out simulations to test 
the equilibrium behavior of the model in the case of the 
isothermal model. This is a crucial test as the deriva- 
tion of the fluctuating forces relies on the most signif- 
icant approximations. The simulations are carried out 
using a periodic Voronoi tesselation described in detail 
elsewhere ]2?| . 



For most practical applications involving complex flu- 
ids, additional interactions and boundary conditions need 
to be specified. These too must be deduced from the mi- 
croscopic dynamics, just as we have done for the inter- 
particle forces. This may be achieved by considering a 
particulate description of the boundary itself and includ- 
ing molecular interactions between the fluid MD parti- 
cles and other objects, such as particles or walls. Ap- 
propriate modifications can then be made on the basis of 
the momentum- flux tensor of Eq. (|l^) , which is generally 
valid. 

Consider for example the case of a colloidal suspension, 
which is shown in Fig. |6| Beginning with the hydrody- 
namic momentum- flux tensor Eq. (|l9|) and Eq. (B9), it is 
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evident that we also need to define an interaction re- 
gion where the DP-colloid forces act: the DP-colloid 
interaction may be obtained in the same form as the 
DP-DP interaction of Eq. (^) by making the replace- 
ment Iki — > Lki, where L k i is the length (or area in 
3D) of the arc segment where the dissipative particle 
meets the colloid (see Fig. ^) and the velocity gradient 
r kl ((Ufci • 6fci)efci + Ufci) is that between the dissipative 
particle and the colloid surface. The latter may be com- 
puted using Ufc and the velocity of the colloid surface 
together with a no-slip boundary condition on this sur- 
face. In Eq. ( |6^ ) the replacement l k i — > Lki must also 
be made. 




FIG. 6. Multiscale modeling of colloidal fluids. As 

usual, the dissipative particles are defined as cells in the 
Voronoi lattice. Note that there are four relevant length scales 
in this problem: the scale of the large, gray colloid particles, 
the two distinct scales of the dissipative particles in between 
and away from the colloids and finally the molecular scale of 
the MD particles. These mediate the mesoscopic interactions 
and are shown as dots on the boundaries between dissipative 
and colloidal particles. 

Although previous DPD simulations of colloidal fluids 
have proved rather successful jlO| at low to intermedi- 
ate solids volume fractions, they break down for dense 
systems whose solids volume fraction exceeds a value of 
about 40% because the existing method is unable to han- 
dle multiple lengthscale phenomena. However, our new 
version of the algorithm provides the freedom to define 
dissipative particle sizes according to the local resolution 
requirements as illustrated in Fig. ^. In order to increase 
the spatial resolution where colloidal particles are within 



close proximity it is necessary and perfectly admissible to 
introduce a higher density of dissipative particles there; 
this ensures that fluid lubrication and hydrodynamic ef- 
fects are properly maintained. After these dissipative 
particles have moved it may be necessary to re-tile the 
DP system; this is easily achieved by distributing the 
mass and momentum of the old dissipative particles on 
the new ones according to their area (or volume in 3D). 
Considerations of space prevent us from discussing this 
problem further in the present paper, but we plan to re- 
port in detail on such dense colloidal particle simulations 
using our method in future publications. We note in 
passing that a wide variety of other complex systems ex- 
ist where modeling and simulation are challenged by the 
presence of several simultaneous length scales, for exam- 
ple in polymeric and amphiphilic fluids, particularly in 
confined geometries such as porous media pfj| . 

B. The low viscosity limit and high Reynolds 
numbers 

In the kinetic theory derived by Marsh, Backx and 
Ernst [15] the viscosity is explicitly shown to have a ki- 
netic contribution tjx — pD/2 where D is the DP self 
diffusion coefficient and p the mass density. The kinetic 
contribution to the viscosity was measured by Masters 
and Warren ]3l|] within the context of an improved the- 
ory. How then can the viscosity rj used in our model 
be decreased to zero while kinetic theory puts the lower 
limit t\k to it? 

To answer this question we must define a physical way 
of decreasing the MD viscosity while keeping other quan- 
tities fixed, or, alternatively rescale the system in a way 
that has the equivalent effect. The latter method is 
preferable as it allows the underlying microscopic sys- 
tem to remain fixed. In order to do this we non- 
dimensionalize the DP momentum equation Eq. (|39|). 

For this purpose we introduce the characteristic equi- 
librium velocity, Uq = \J ksT/M, the characteristic dis- 
tance ro as the typical DP size. Then the characteristic 
time t' = ro/Uo follows. 

Neglecting gravity for the time being Eq. ([39]) takes 
the form 

^ = - E *« (if eH + Fie" (u ' fci + (u ^ ■ ekl)ea) ) 

I \ " l 'kl L 'kl Pk + Pi - it! U fc + U i i V* ■&/ / fi o\ 

+ Z^^v 2 — kl 2 + 2^ b ki, (o») 

i kl i 

where = P k /(MU ), p' H = p kl rl/(MU$), M = pr 2 Q 
in 2d, the Reynolds number Re = Uor p/rj and F' fe; = 
(r /MU$)F kl where F kl is given by Eqs. @ and @. 
A small calculations then shows that if F' kl is related to 
u)' kl and t' like F k i related to uj k i and t, then 
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.'2 



1 k B T 
Re MXJl 



1 

Re 



(69) 



where we have neglected dimensionless geometric prefac- 
tors like Iki/fkl and used the fact that the ratio of the 
thermal to kinetic energy by definition of Uq is one. 

The above results imply that when the DPD system is 
measured in non-dimcnsionalizcd units everything is de- 
termined by the value of the mesoscopic Reynolds num- 
ber Re. There is thus no observable difference in this 
system between increasing r$ and decreasing rj. 

Returning to dimensional units again the DP diffusiv- 
i ty may be obtained from the Stokes-Einstein relation 
[B2| as 



D 



k B T 
ar or/ 



(70) 



where a is some geometric factor (a — 6tt for a sphere) 
and all quantities on the right hand side except rg refer 
directly to the underlying MD. As we are keeping the 
MD system fixed and increasing Re by increasing ro, it 
is seen that D and hence rjx vanish in the process. 

We note in passing that if D is written in terms of the 
mean free path A: D = X^JksT / (pr§) and this result is 
compared with Eq. ( |70| ) we get A' = X/r ~ l/r in 2d, 
i.e. the mean free path, measured in units of the particle 
size decreases as the inverse particle size. This is consis- 
tent with the decay of t\k- The above argument shows 
that decreasing r/ is equivalent to keeping the microscopic 
MD system fixed while increasing the DP size, in which 
case the mean free path effects on viscosity is decreased 
to zero as the DP size is increased to infinity. It is in this 
limit that high Re values may be achieved. 

Note that in this limit the thermal forces Fm ~ Re -1 / 2 
will vanish, and that we arc effectively left with a macro- 
scopic, fluctuationless description. This is no problem 
when using the present Voronoi construction. However, 
the effectively spherical particles of conventional DPD 
will freeze into a colloidal crystal, i.e. into a lattice config- 
uration [8,9] in this limit. Also while conventional DPD 
has usually required calibration simulations to determine 
the viscosity, due to discrepancies between theory and 
measurements, the viscosity in this new form of DPD is 
simply an input parameter. However, there may still be 
discrepancies due to the approximations made in going 
from MD to DPD. These approximations include the lin- 
earization of the inter-DP velocity fields, the Markovian 
assumption in the force correlations and the neglect of a 
DP angular momentum variable. 

None of the conclusions from the above arguments 
would change if we had worked in three dimensions in 
stead of two. 



VII. CONCLUSIONS 

We have introduced a systematic procedure for de- 
riving the mesoscopic modeling and simulation method 
known as dissipative particle dynamics from the under- 
lying description in terms of molecular dynamics. 



MD 



coarse graining 



Time symmetric 
coarsed grained equations 
of motion 



| averaging over MD configurations 



DPD1 

| constitutive relations and Markovian assumption 



DPD2: Langevin equations 



Fokker-Planck equations 

1 equilibrium statistical mechanics 




Fluctuating 
hydrodynamics 



Fluctuation dissipation co = co(r| , t ) 
DPD complete 



FIG. 7. Outline of the derivation of dissipative particle dy- 
namics from molecular dynamics as presented in the present 
paper. The MD viscosity is denoted by r\ and ui is the ampli- 
tude of the fluctuating force F as defined in Eq. ( fii| ) 

Figure Q illustrates the structure of the theoretical de- 
velopment of DPD equations from MD as presented in 
this paper. The initial coarse graining leads to equa- 
tions of essentially the same structure as the final DPD 
equations. However, they are still invariant under time- 
reversal. The label DPD1 refers to Eqs. ©, (|J) and 
(fiof), whereas the DPD2 equations have been supple- 
mented with specific constitutive relations both for the 
non-equilibrium fluxes (momentum and heat) and an 
equilibrium description of the thermodynamics. These 
equations are Eqs. ( |39| ) and ( |4g| ) along with Eqs. (|6^). 
The development we have made which is shown in Fig. |t] 
does not claim to derive the irreversible DPD equations 
from the reversible ones of molecular dynamics in a rig- 
orous manner, although it does illustrate where the tran- 
sition takes place with the introduction of molecular av- 
erages. The kinetic equations of this new DPD satisfy 
an iJ-theorem, guaranteeing an irreversible approach to 
the equilibrium state. Note that in passing to the time- 
asymmetric description by the introduction of the av- 



eraged description of Eq. (36), a time asymmetric non- 
equilibrium ensemble is required j23| . 

This is the first time that any of the various exist- 
ing mesoscale methods have been put on a firm 'bot- 
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torn up' theoretical foundation, a development which 
brings with it numerous new insights as well as prac- 
tical advantages. One of the main virtues of this proce- 
dure is the capability it provides to choose one or more 
coarse-graining lengthscales to suit the particular model- 
ing problem at hand. The relative scale between molecu- 
lar dynamics and the chosen dissipative particle dynam- 
ics, which may be defined as the ratio of their number 
densities pdpd/pmd, is a free parameter within the the- 
ory. Indeed, this rescaling may be viewed as a renormal- 
isation group procedure under which the fluid viscosity 
remains constant: since the conservation laws hold ex- 
actly at every level of coarse graining, the result of doing 
two rescalings, say from MD to DPDa and from DPDa 
to DPD/3, is the same as doing just one with a larger 
ratio, i.e. pDPDfi/pMD = (pdpd/j/ PdpdqX/Odpdo/ Pmd)- 

The present coarse graining scheme is not limited to 
hydrodynamics. It could in principle be used to rescale 
the local description of any quantity of interest. However, 
only for locally conserved quantities will the DP particle 
interactions take the form of surface terms as here, and 
so it is unlikely that the scheme will produce a useful 
description of non-conserved quantities. 

In this context, we note that the bottom-up approach 
to fluid mechanics presented here may throw new light 
on aspects of the problem of homogeneous and inhomo- 
geneous turbulence. Top-down multiscale methods and, 
to a more limited extent, ideas taken from renormali- 
sation group theory have been applied quite widely in 
recent years to provide insight into the nature of turbu- 
lence |53|,|4j ; one might expect an alternative perspective 
to emerge from a fluid dynamical theory originating at 
the microscopic level, in which the central relationship 
between conservative and dissipative processes is speci- 
fied in a more fundamental manner. From a practical 
point of view it is noted that, since the DPD viscosity is 
the same as the viscosity emergent from the underlying 
MD level, it may be treated as a free parameter in the 
DPD model, and thus high Reynolds numbers may be 
reached. In the 77 — > limit the model thus represents 
a potential tool for hydrodynamic simulations of turbu- 
lence. However, we have not investigated the potential 
numerical complications of this limit. 

The dissipative particle dynamics which we have de- 
rived is formally similar to the conventional version, in- 
corporating as it does conservative, dissipative and fluc- 
tuating forces. The interactions are pairwise, and con- 
serve mass and momentum as well as energy. However, 
now all these forces have been derived from the under- 
lying molecular dynamics. The conservative and dissi- 
pative forces arise directly from the hydrodynamic de- 
scription of the molecular dynamics and the properties of 
the fluctuating forces are determined via a fluctuation- 
dissipation relation. 

The simple hydrodynamic description of the molecules 
chosen here is not a necessary requirement. Other choices 



for the average of the general momentum and energy flux 
tensors Eqs. ( p6| ) and ( |l9| ) may be made and we hope 
these will be explored in future work. More significant 
is the fact that our analysis permits the introduction of 
specific physicochemical interactions at the mesoscopic 
level, together with a well-defined scale for this meso- 
scopic description. 

While the Gaussian basis we used for the sampling 
functions is an arbitrary albeit convenient choice, the 
Voronoi geometry itself emerged naturally from the re- 
quirement that all the MD particles be fully accounted 
for. Well defined procedures already exist in the litera- 
ture for the computation of Voronoi tesselations |3j| and 
so algorithms based on our model are not computation- 
ally difficult to implement. Nevertheless, it should be ap- 
preciated that the Voronoi construction represents a sig- 
nificant computational overhead. This overhead is of or- 
der iVlog N, a factor logiV larger than the most efficient 
multipolc methods in principle available for handling the 
particle interactions in molecular dynamics. However, 
the prefactors are likely to be much larger in the particle 
interaction case. 

Finally we note the formal similarity of the present par- 
ticulate description to existing continuum fluid dynam- 
ics methods incorporating adaptive meshes, which start 
out from a top-down or macroscopic description. These 
top-down approaches include in particular smoothed par- 
ticle hydrodynamics |l^] and finite-element simulations. 
In these descriptions too the computational method is 
based on tracing the motion of elements of the fluid on 
the basis of the forces acting between them |36|. However, 
while such top-down computational strategies depend on 
a macroscopic and purely phenomenological fluid descrip- 
tion, the present approach rests on a molecular basis. 
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